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We focus on a multidimensional field with uncorrelated spectrum, and study the quality of the 
reconstructed signal when the field samples are irregularly spaced and affected by independent and 
identically distributed noise. More specifically, we apply linear reconstruction techniques and take the 
mean square error (MSE) of the field estimate as a metric to evaluate the signal reconstruction quality. 
Q , We find that the MSE analysis could be carried out by using the closed-form expression of the eigenvalue 

distribution of the matrix representing the sampling system. Unfortunately, such distribution is still 
J> ■ unknown. Thus, we first derive a closed-form expression of the distribution moments, and we find 

in 1 

ly-j ' that the eigenvalue distribution tends to the Marcenko-Pastur distribution as the field dimension goes to 
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infinity. Finally, by using our approach, we derive a tight approximation to the MSE of the reconstructed 



■^J- . field. 

o 

OO 

O 



I. Introduction 



We address the important issue of reconstructing a multidimensional signal from a collection of samples 
that are noisy and not uniformly spaced. As a case study, we consider a wireless sensor network for 
environmental monitoring, where the nodes sensing the physical phenomenon (hereinafter also called 
field) are randomly deployed over the area under observation. The sensors sample a d-dimensional 
spatially finite physical field, where d may take into account spatial dimensions as well as the temporal 
dimension. Examples of such fields are pressure or temperature, on a 4-dimension domain, i.e., three 
spatial coordinates plus the time dimension. A spatially finite physical field is not bandlimited, however it 
admits an infinite Fourier series expansion. Here, we consider a finite approximation of the physical field 
obtained by truncating such series, assuming that the contribution of the truncated terms is negligible. 

This work was supported by MIUR through the MEADOW project 
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In our case study, we assume that the measured samples are transferred from the sensors to a common 
data-collecting unit, the so-called sink node, which is in charge of reconstructing the field. Distributed 
systems, where in-network processing, is performed are out of the scope of this work. We do not deal 
with issues related to information transport and, thus, we assume that all samples are correctly received 
at the sink node. The field samples, however, are corrupted by additive i.i.d. noise, due to quantization, 
round-off errors or quality of the sensing device. Furthermore, the sampling points are known at the sink 
node, because (i) either sensors are located at pre-defined positions or their position can be estimated 
through a localization technique [1], and (ii) the sampling time is either periodic or included in the 
information sent to the sink. 

Several efficient and fast algorithms have been proposed to numerically reconstruct or approximate 
a signal in such setting, which amount to the solution of a linear system (see [2], [3] and references 
therein). A widely used technique consists in processing the sensors' measures by means of a linear filter, 
which is a function of the system parameters known at the sink. We observe that the following two major 
factors affect the linear reconstruction: 

(i) the given machine precision, which may prevent the reconstruction algorithm from performing 
correctly and may lead to a non-negligible probability of reconstruction failure [3], 

(ii) the noise level affecting the sensors' measurements. 

In the latter case, a measure of reconstruction accuracy is given by the mean square error (MSE) of 
the field estimate. In [4], [5], we have found that these issues could be studied by using the eigenvalue 
distribution of the reconstruction matrix; however, obtaining such a distribution is still an open problem. 
In this work, we first extend the system model and the problem formulation presented in [5] to the case 
of multidimensional fields (Section ITITt - Then, we derive a closed-form expression of the moments of 
the eigenvalue distribution, through asymptotic analysis (Section IV]). By using the moments expressions, 
we prove that the eigenvalue distribution of the matrix representing the sampling system tends to the 
Marcenko-Pastur distribution [6] as the field dimension d —> oo (Section [VTt . 

We apply our results to the study of the MSE of the field estimate, when the sensors measurements 
are noisy and the reconstruction at the sink is performed through linear filtering. 

We generalize the MSE expressions to the multidimensional case (with finite d), and we show that, 
by using the Marcenko-Pastur distribution instead of the actual eigenvalue distribution, we obtain an 
approximation to the MSE of the reconstructed field which is very tight for d > 3 (Section I VIII ). 

Before providing a detailed description of our analysis, in the next section we discuss some related 
studies and highlight our main contributions with respect to previous work. 
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II. Related work and main contributions 

Relevant to our work is the literature on spectral analysis, where, however, several studies deal with 
regularly sampled signals (e.g., [7] and references therein). An excellent guide to irregular sampling is 
[8], which presents a large number of techniques, algorithms, and applications. Reconstruction techniques 
for irregularly or randomly sampled signals can be found in [3], [9], [10], just to name few. In particular, 
Feichtinger and Grochenig in [10] provide an error analysis of an iterative reconstruction algorithm taking 
into account round-off errors, jitters, truncation errors and aliasing. From the theoretical point of view, 
irregular sampling has been studied in [3], [9]-[14] and references therein. 

In the context of sensor networks, efficient techniques for spatial sampling are proposed in [15], [16]. 
In particular, in [16], an adaptive sampling is described, which allows the central data-collector to vary 
the number of active sensors, i.e., samples, according to the desired resolution level. Data acquisition is 
also studied in [17], where the authors consider a unidimensional field, uniformly sampled at the Nyquist 
frequency by low precision sensors. The authors show that the number of samples can be traded-off with 
the precision of sensors. The problem of the reconstruction of a bandlimited signal from an irregular set 
of samples at unknown locations is addressed in [18]. There, different solution methods are proposed, 
and the conditions for which there exist multiple solutions or a unique solution are discussed. Differently 
from [18], we assume that the sink can either acquire or estimate the sensor locations and that sensors 
are randomly deployed. 

The field reconstruction at the sink node with spatial and temporal correlation among sensor measures 
is studied, for instance, in [19]— [23]. Other interesting studies can be found in [24], [25], which address 
the perturbations of regular sampling in shift-invariant spaces [24] and the reconstruction of irregularly 
sampled images in presence of measurement noise [25]. 

We point out that our main contribution with respect to previous work on signal sampling and 
reconstruction is the probabilistic approach we adopt to analyze the quality level of a signal reconstructed 
from a set of irregular, noisy samples. Our analysis, however, applies to sampling systems where the 
field reconstruction is performed in a centralized manner. Finally, we highlight that our previous work [5] 
assumes that sensors are uniformly distributed over the spatial observation interval and may be displaced 
around a known average location. The effects of noisy measures and jittered positions are analyzed when 
linear reconstruction techniques are employed. However, only the unidimensional case is studied and 
semi-analytical derivations of the MSE of the reconstructed field are obtained. In [26], instead, sensors 
are assumed to be fixed, and the objective is to evaluate the performance of a linear reconstruction 
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technique in the presence of quasi-equally spaced sensor layouts. 
A. Main results 

The goal of this work is to provide an analytical study on the reconstruction quality of a multidi- 
mensional physical field, with uncorrected spectrum. The field samples are (i) irregularly spaced, since 
they are gathered by a randomly deployed sensor network and (ii) affected by i.i.d. noise. The sink 
node receives the field samples and runs the reconstruction algorithm in a centralized manner. Our major 
contributions with respect to previous work are as follows. 

1 . Given a (i-dimensional problem formulation, we obtain analytical expressions for the moments of the 
eigenvalue distribution of the reconstruction matrix. Using the expressions of the moments, we show 
that the eigenvalue distribution tends to the Marcenko-Pastur distribution [6] as the field dimension 
d — > oo. 

2. We apply our results to the study of the quality of a reconstructed field and derive a tight approxi- 
mation to the MSE of the estimated field. 

III. Preliminaries 

We first present the multidimensional formulation of our reconstruction problem. Then, we give some 
background on linear reconstruction techniques and generalize to the multidimensional case some results 
previously obtained in the unidimensional case [5]. Finally, we highlight the main steps followed in our 
study. 

Notation: Lower case bold letters denote column vectors, while upper case bold letters denote matrices. 
We denote the (h, k)-th entry of the matrix X by (X)^, the transpose of x by x T , and the conjugate 
transpose of x by xt The identity matrix is denoted by I. Finally, E[x] is the average of x and subscripts 
to the average operator specify the variable with respect to which the average is taken. 

A. Irregular sampling of multidimensional signals 

Let us consider a <i-dimensional, spatially-finite physical field (d > 1), where r sensors are located 
in the hypercube H = {x | x G [0, l) d } and measure the value of the field. We assume that the sensor 
sampling points are known. At first, we consider that they are deterministic, then we will assume that 
they are i.i.d. random variables uniformly distributed in the hypercube H. 
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When observed over a finite region, a d-dimensional physical field s(x) with finite energy E s admits 
an infinite d-dimensional Fourier series expansion with coefficients an, such that E s = E^°° t d =-oo I^P 
where £ = [£i, . . . , £J\ is a vector of integers and l m , m = 1, . . . , d represents the index of the expansion 
along the m-th dimension. We truncate the expansion to 2M + 1 terms per dimension where M is such 
that ^7i M ~£d=-oo l^l 2 + E^°° i d =M+i \^t\ 2 ^ E s Therefore, one can think of M as the approximate 
one-sided bandwidth (per dimension) of the field, which can be approximated over the finite region H 
as 

S ( X ) = (2M + iy d ' 2 ^(^ 2 " T£ (i) 

I 

where the term (2M + l)~ d / 2 is a normalization factor and Yle. represents a d-dimensional sum over the 
vector £, with i m = —M, ... ,M. Also, a u ^ = di and the function 

d 

u(£) = Y,m+i) m - i £ m , 

m=l 

_(2M-kl)_-i ^ v ^ ^ _^_ (2M+i) -l ma p S ^ vector £ on t a sca i ar index. Note that, while dg has a 
vectorial index, a v tg\ has a scalar index and it has been introduced to simplify the notation. As an example, 
for d = 2 and M = 1, we have z/(£) = 34+£> and a(xi,x 2 ) = 5 E< 1= -i E/ 2 =-i o % +/ 1 e i2,r ( a!l<1+a!a<! '). 

Let A? = {xi, . . . , x r }, with x g = [a^i, . . . , x qt d] T £ H, q = 1, . . . , r, be the set of sampling points, 
and s = [si, . . . , s r ] T , s q = s(x g ), the values of the corresponding field samples. Following [3], we write 
the vector of field values s as a function of the spectrum: 

s = G+a (2) 

where a is a vector of size (2M + l) d , whose u(£)-th entry is given by a u ^, and G^ is the (2M + l) d x r 
matrix: 

{G d ) u{tu = {2M + l)- d '\-^ 1 (3) 

In general, the entries of a can be correlated with covariance matrix E[aa^] = a^C a , and Tr{C a } = 
(2M + l) d . In the following, we restrict our analysis to the class of fields characterized by E[aa^] = a^I. 
If the sensor measurements, p = [pi, . . . ,p r } T , are noisy, then the relation between sensors' samples and 
field spectrum can be written as: 

P = s + n = G^a + n (4) 

where the noise is represented by the r-size, zero-mean random vector n, with covariance matrix E[nnt] = 
a%L r . We define the signal-to-noise ratio on the measure as: SNR m ^ al/a 2 n = 1/a. 
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B. Sampling rate 

Following [5], we introduce the parameter f3 denned as: 

0= M±iH (5) 

r 

This parameter represents the ratio between the number of harmonics used for the field reconstruction 
and the number of sensors sampling the field. In the following, we consider f3 G [0, 1). Notice that for 
fixed and M, the number r of samples exponentially increases with d. 

C. Previous results on reconstruction quality 

Given an estimate a of the field spectrum a, the reconstructed signal is: 

S(x) = (2M + iy d ' 2 £ ^(V^ (6) 
e. 

As reconstruction performance metric, we consider the MSE of the field estimate, which, for any given 
set of sampling points X, is defined as: 

E ||a-a|| 2 

MSE, = E I |S(x) - 3 (x)P dx = (7, 

where the average is taken with respect to the subscripted random vectors. Note that ([TJ) still assumes 
that the sampling points are deterministic; this assumption will be removed later in the paper. 

For linear models such as (0]), several estimation techniques based on linear filtering are available in 
the literature [27]. We employ a filter B such that the estimate of the field spectrum is given by the 
linear operation 

a = B f p (8) 

where B is an r x (2M + l) d matrix. In particular, we consider the linear filter providing the best 
performance in terms of MSE, i.e., the linear minimum MSE (LMMSE) filtefl [27]: 



B = Gt(Rd + aiy 1 (9) 



where K d = G d G\. 



From now on, we carry out our analysis under the assumption that the elements of the set X are 
independent random vectors, with i.i.d. entries, uniformly distributed in the hypercube 7i. 

In [5], we have shown that a simple and effective tool to evaluate the performance of large finite 
systems is asymptotic analysis. We computed the MSE by letting the field number of harmonics and the 

'Notice that when the covariance matrix of a is known, the LMMSE filter generalizes to (G^CaGj + oJ) -1 G^C a . 
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number of samples grow to infinity, while their ratio (3 = (2M + l) d /r is kept constant. We observed 
the validity of asymptotic analysis results, even for small values of M and r. Similarly, here we consider 
as performance metric the asymptotic average MSE, normalized to <7„: 

1 



MSE r 



M,r- 



lim 



■ e[mse 



(10) 



where below the limit denotes the ratio which is kept constant. In (fTOl . the average is over all possible 
realizations of the set X. Using ©-(O and the above definitions, in Appendix |E] we show that, in the 
case of the LMMSE: 



MSE r 



E 

*d,0 



a/3 



en) 



\ dij 3 + a(3_ 

where is a random variable with probability density function (pdf) fd,p{x), distributed as the 
asymptotic eigenvalues of = (3Hd = (3GdG\. The subscripts d and f3 of A indicate that the distribution 
of the asymptotic eigenvalues of depends on both the field dimension d and the parameter (3. 

The matrix plays an important role in our analysis; in the following, we introduce some of its 
properties. In the unidimensional case (d = 1), Ti is a (2M + 1) x (2M + 1) Hermitian Toeplitz matrix 
given by 

( to h ■ ■ ■ t 2 M ^ 

t-i to • • • t2M~X 

Tj = 

\ t-2M t-2M+l - - - *0 / 

where (Tx)^ = t t _ t . = \ exp(-j27r(^ - l')x q ), £,£' = -M,...,M. For d > 2, T d can be 

defined recursively as a (2M + l) d x (2M + l) d Hermitian Block Toeplitz matrix with non Hermitian 
Toeplitz blocks: 

/ B Bi • • • B 2 Af ^ 

B_i B • • • B 2 m-i 



B 



-2M+1 



B, 



where 



-d)v(£),is(£') 



1 r 

q=l 



-}2n(e.-£')x. q 



(12) 



and £,£' £ [-M, ...,M] d . That is, the matrix T d is composed of (2M + l) 2 blocks Bj of size 
(2M + x (2M + each including (2M + l) 2 blocks of size (2M + l) d ~ 2 x (2M + l) d ~ 2 , 
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and so on. The smallest blocks have size (2M + 1) x (2M + 1); they have the same structure as matrix 
Ti in the unidimensional case, however only those on the main diagonal have a Hermitian structure. 
Proof of this is given in [28] for d = 2; the extension to the d-dimensional case is straightforward. 

IV. Estimation Error Calculation Method 
The analysis detailed in the next sections consists of the following main steps. 

(i) As a practical case, we consider the asymptotic expression of the LMMSE in (fTTT > and notice that 
an analytical evaluation of the asymptotic LMMSE could be obtained by exploiting the closed- 
form expression of the eigenvalue distribution, fd,p(x), of the reconstruction matrix. However, such 
expression is still unknown. Hence, as a first step we derive a closed form expression of the moments 
of Arf p, for any d and (3, and provide an algorithm to compute them. 

( ii) We show that the value of the moments of the eigenvalue distribution decreases as the field dimension 
d increases. 

(Hi) We prove that, as d — > oo, the expression of the eigenvalue distribution tends to the Marcenko-Pastur 
distribution. 

(iv) By using the Marcenko-Pastur distribution, we are able to obtain a tight approximation for the 
LMMSE of the reconstructed field, which holds for any finite value of d. 

V. Closed Form Expression of the moments of the asymptotic eigenvalue pdf 

Ideally, we would like to obtain the analytical expression of the distribution fd,p(x) of the asymptotic 
eigenvalue of T^, for a given (3. Unfortunately, such a calculation seems to be prohibitive and is still an 
open problem. Therefore, as a first step, we compute the closed form expression of the moments E[A^ J 
of \d,p, for any positive integer p. 

In the limit for M and r growing to infinity with constant j3, the expression of E[A^ J can be easily 
obtained from the powers of as in [29], [30], 

E[A?J= lim — i — -rTr E [T^l (13) 

L W\ M.r-^+oo (2M + 1 ) d X L di 

P 

In Section IV-AI we show that E[A^ ^] is a polynomial in 0, of degree p — 1 (see d24b); the remaining 
subsections describe how to compute this polynomial. 
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A. Partitions 



Using CED, the term TrE[T^] in C[3]> can be written as 



x 



TrEM 



x 



E 

,Y 



5j( T d)i'(ii),w(^) 



E 

A' 



^ ( T rf)^(^l),^(£2) • • • ( T ^ 



1 



(14) 



f 1 T 



[< + l] 



LLe£d 

E E | 

qefi LeCd 

where Q = {q | q = . . . , q p ]}, q { = 1, . . . , r C d = {L | L = . . . , = 
ii_ = —M, ... M and 

i + 1 1 < % < p 
1 i = p 

In (fl4l . the average is performed over the random set of positions X = {xi, . . . , x r }, with independent 
and uniformly distributed elements. To obtain a closed-form expression of the distribution moments, we 
rewrite ([141 by using set partitioning. 

Let V = {1, . . . ,p} be the set of integers from 1 to p. We observe that any given vector q G Q partitions 
the set V into 1 < A;(q) < p disjoint non-empty subsets 'Pi(q), . . . ,7- > fe(q), where Vj, j = 1, . . . , k(q), 
is the set of indices of the entries of q taking the same value jj. That is, 



Vj(q) = {ieV | qi = lj} 



(15) 



and A;(q) is the number of distinct values jj taken by the entries of vector q. Subsets Vj have the 
following properties 

fc &V-(q)=P, ^(q)n^(q)=0 
J=l 

for, j 7^ j'. Also, we point out that, since r is the number of values that the entries % can take, there 
exist r\/(r — k(q))\ vectors q E Q generating a given partition of V made of k(q) subsets. In order to 
clarify the above concepts, we provide an example below. 
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Example 1: Let p = 6, then V = {1,2,3,4,5,6}. Also, let q = [4,9,5,5,4,3]. Since the 
distinct values in q are jj = 3, 4, 5, 9, we have /c(q) = 4. It follows that V is partitioned into 
the following subsets: 

7>i(q) = {l,5} (91 = 95 = 4) 
7> 2 (q) = {2} (92 = 9) 
P 3 (q) = {3,4} ( 93 = ?4 = 5) 
^(q) = {6} (96 = 3) 
Hence, the partition of V induced by q is {{1, 5}, {2}, {3,4}, {6}}. 



Next, we introduce an effective method to represent partitions of a set V, by building a tree of depth 
p, as in Figure Q] Such a representation will allow to simplify the notation in the following analysis. To 
build the tree of depth p, we proceed as follows. Each node of the tree is assigned with a label from 
the set V = {1, . . . , p}, starting from the root which is labeled by 1. Each node generates m + 1 leaves, 
labeled in increasing order from 1 to m + 1, where m is the largest label on the path from the root to 
such node. Note that, at level p, any value in {1, . . . ,p} is used to label the leaves at least once. 

Then, given a tree of depth p, we define u> = [uj\, . . . , u> p ] as a path of length p from the tree's root 
to a leaf. We observe that a vector q can be represented as a path cj in the tree of depth p. This is done 
by assigning a label (j = 1, . . . ,p) in increasing order to every distinct value of q; the vector collecting 
the labels is the path u corresponding to the given q. We have that the path cj, corresponding to a given 
q, defines in the tree of depth p the same partition of the set V as the one induced by q. Indeed, given 
a partition of V, the subset Vj defined in ( fT5T ) can be rewritten as 

P j (u;) = {ier | tJi = j}, (16) 

i.e., as the set of integers corresponding to the depths of the j-th label in the path. As a last remark, 
consider the number of distinct values in q (i.e., the number of distinct labels in u) to be equal to fc(q), 
and recall that the number of all possible values taken by the p elements of q is equal to r. It follows 
that r!/(r — k(q))l different q's yield the same vector cj. This is in agreement with the fact that, given 
V, there are r!/(r — £;(q))! different q's generating the same partition consisting of k(q) subsets. Again, 
for the sake of clarity, we give an example. 
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Fig. 1. Partitions tree of depth p — 4. The path uj = [1, 2, 1, 1] employed in Example 2 is highlighted by using dashed lines 



Example 2: Let us consider p = 4 and q = [4,9,4,4]. The vector q can be represented in 
the tree of depth 4 as the path u> = [1, 2, 1, 1] (i.e., the path highlighted with dashed lines in 
Figure Q]). In q there are two distinct values (namely, 4 and 9), or, equivalently, in the path 
<jj there are two labels (namely, 1 and 2); then k(q) = k(u) = 2. Label 1 appears in u at 
depths 1,3, and 4 (V\ = {1, 3, 4}), while label 2 appears at depth 2 (V2 = {2}). The partition 
of V = {1, 2, 3, 4} induced by q or, equivalently, by u is therefore: {{1, 3, 4}, {2}}. 

From the discussion above, it should be clear that considering a partition of V is equivalent to 
considering a path u in a tree of depth p. Hence, in the following analysis, we will refer to a partition 
through its corresponding path u. 

We now exploit set partitioning to rewrite (fT4l . Since the random vectors x g are independent, given 
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q, the average operator in (fl4l ) factorizes into fc(q) terms, i.e., 

a -j»rElLi*? 4 (A-^+i]) 



E 

X 



fc(q) 



(17) 



Each term depends on a single random vector x T . Moreover, since the entries of x 7j are independent 
random variables uniformly distributed in [0, 1), we have: 



E 



-j27TxT £ ie p, (q) 



1 [ e [ e - j27r ^- c ^i = n 5 ( 



m=l Tj 



(18) 



m=l 



where x~ /jj7n and 4,m are the m-th entries of x T and t^, respectively, where the function 5(-) is the 
Kronecker's delta, and where Cj m = YlieV M ^*> m ~ %+i],m- By substituting (fTTb and ( fT8l ) in (fl4l ) and 
by expanding the summation ^Le£d' we obtain 

., fc(q) d 

v, e e e n n *<<*-> 



TrE[T>] 



qeQ £iG£i £ d e£i i=l m=l 



1 d 

e n 



qeQ m=l 



1 E 



qeQ 



fc(q) 

e n*( 

m e£i j=i 



fc(q) 

eii%) 

£e£i j=i 



(19) 



where Cj = YlieV-M ^ ~ F° r an y gi ven q £ Q, the expression 

fc(q) 

cm(q) = e n %) ( 2 °) 

£e£i i=i 

is a polynomial in M, since it represents the number of points with integer coordinates contained in the 
hypercube [— M, . . . , M] p and satisfying the k(q) constraints: 

c i= E ^-%+i]=0 (21) 

for j = 1, ... , fc(q). In Appendix lAl we show that one of these constraints is always redundant and that the 
number of linearly independent constraints is exactly equal to k(q) — 1. As a consequence, the polynomial 
Ciw(q) has degree p — fe(q) + 1, and, for large values of M, we have Cm(q) = 0((2M + l) p ~ fc ( q ) +1 ). 
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Now, using ( fl9l ) and (1201 ). the limit in (TT3T > is given by 

- „«5U S <22) 

Equation (|22l) can be further simplified by considering that there exist r!/(r — fc(q))! vectors q 6 Q 
generating a given partition of V made of fc(q) subsets, or, equivalently, a path oj of length p with fc(q) 
distinct labels. 

Let O p be the set of vectors u>, each corresponding to a distinct partition of T 3 . Also, let us write k(q) 
and CA/(q) as functions of the partition induced by q, i.e., as functions of u>. Then, we obtain: 

= lim V Cm (faj) d r! 

where 

• the notation Ylq=>cj represents the sum over all vectors q generating a given path u>, 

• the equality (a) holds because the number of vectors q generating a given u is r!/(r — A;(u;))!. 
Note that, for large r, r\/(i — &(>>))! = r*( w ' + 0{r k ^~ l ). Also, since Cm (^) is a polynomial in M 
of degree p - k(u) + 1, for large values of M we have: Cm O) = f(u;)(2M + + 0((2M + 
l)P- fc H), where v(u) is the coefficient of degree p — k(u) + 1 of Cm(<^)- Therefore, taking the limit, 
we obtain: 

e[au = E vWf-"™ = E cm) 

wefip fc=i ue$] Pi t 

where Q, p % C P is the subset of O p containing paths with k(u>) distinct labels, and 

„(„) = lim 6"^) ^ (25) 

v ; M-»+oo (2M + l)P-fc(«)+i 

Note that the coefficient v(u>) represents the volume of the convex poly tope described by the constraints 

in (l2Tb . when the variables ti are considered as real and limited to a j>-dimensional hypercube of volume 

1. As a consequence, we have: < v{u) < 1. 

Equation (l24l ) provides a closed-form expression of the moment E[A^ J, as a polynomial in (3 of degree 

p — 1. Again, for the sake of clarity, we give an example below. 
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Example 3: Let p = 6 and q = [4, 9, 5, 5, 4, 3]. We have 
of V is {{1, 5}, {2}, {3, 4}, {6}}. Then, the set of k(u) 



= [1, 2, 3, 3, 1, 4], and the partition 
4 constraints (as in (1211 ) are given 



by: 



e 1 + e 5 = i 2 + e 6 



£■ 



£ 3 + h = h + £ 5 



6 



The last equation is redundant since can be obtained from the first three constraints. 



Simplifying, we obtain l\ = 1%, and li = £3 = £5. Since each variable ti ranges from — M to 
M, the number of integer solutions satisfying the constraints is exactly Ca/(u>) = (2M + l) 3 , 
and then v(u>) = lim M ^ +00 (2M+ ^l"li u)+ i = 1. 



Next, in order to compute E[A^ J, we need: 



• to enumerate the partitions, i.e., the vectors cj G Pl fc, for each k = 1, . . . ,p (see Section IV-BI ); 



• to compute the coefficients v(u), for any u> G and k = 1, . . . ,p (see Section |V-C1 ). 

B. Partitions enumeration 

We notice that fi p represents the set of partitions of a p-element set, thus it has cardinality |fi p | = B{p), 
where B(p) is the p-th Bell number or exponential number [31]. Furthermore, the subset fL, C fi p has 
cardinality fc), which is a Stirling number of the second kind [32] given by: 



C. Computation of the coefficients v(u) 

The last step required for the computation of E[A^] is the evaluation of the coefficients v{u), for 
every u £ O p . We have the following Lemma: 

Lemma 5.1: For any u> G fi p (or, equivalently, any partition of V) and any arbitrary integer n, with 
n = 1, . . . , the coefficient u(u;) in (|25T ) is given by: 




with b(p) = ELi s (p^)- 




(26) 



where 
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• yn = [yi, ■ ■ .,y n -i,y n +i, ■ ■ - ,?//cm] t . 

. , , sin(vry) 

. sinc(y) = . 

Try 

Proof: The proof can be found in Appendix |B] ■ 

D. A practical method for the moments computation 

Equation d26l ) in Section IV-CI shows that the computation of E[A^] requires the evaluation of B{p) 
integrals. However, B{p) is very large even for small p, e.g., _B(10) = 115975 and B(20) « 5 • 10 13 . 

The computational complexity can be reduced by recursively applying the simplification rules defined 
in the following Lemma: 

Lemma 5.2: Let 

• u> = [ux, . . . , top] be the path in a tree of depth p, corresponding to the partition of V into k subsets; 

• Vi, . . . ,Vk be the subsets of V = {1, . . . ,p} defined as in (fl6l ); 

• u' be the path obtained from cj by removing cjj. 
We have the following rules: 

1) if Vj has cardinality 1 (i.e., Vj is a singleton) or 

2) if Vj contains adjacencies (in the circular sense), i.e., both i and [i + 1] € Vj, 
then v(u>) = v(u') 

Proof: The proof is a direct consequence of Lemma 15.11 and can be found in Appendix O ■ 
Table U shows two examples of how the rules described in Lemma I5T21 can be applied. Example 1 in the 
Table assumes p = 6 and u> = [1, 2, 3, 2, 2, 1]. At step 1, we note that the third element (i = 3) of u is a 
singleton, then, by applying rule 1, we can remove it from the path. At step 2, we find that in u there 
are some adjacencies, hence we apply twice rule 2 (steps 2 and 3). At step 4, the second element of u 
is a singleton, and we remove it by applying rule 1. Eventually, at step 7, the path lj is empty (i.e., has 
size p = 0) and, thus, the corresponding coefficient is v(u) = 1 (EfA^] = E[l] = 1). 

Example 2 in the Table assumes p = 6 and u> = [1, 2, 3, 1, 2, 1]. After removing a singleton (step 1) 
and an adjacency (step 2), the remaining path cannot be further reduced. Then, to compute the coefficient 
v{u), we need to apply directly Lemma IBTTI on the path u> = [1, 2, 1, 2]. We obtain: 

/+oo 
sinc(yi - y 2 ) sinc(y 2 - Vi) 
-oo 

•sinc(yi - y 2 ) sinc(y 2 - Vi) \y 2 =o dyi 

/+oo 2 
sinc(yi) 4 dyi = - 
-oo 
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TABLE I 

Example of complexity reduction using the rules described in Lemma |5T21 





Example 1 


Example 2 


Step 


u 


Rule 


i 




Rule 


i 


1 


[1,2,3,2,2,1] 


1 


3 


[1,2,3,1,2,1] 


1 


3 


2 


[1,2,2,2,1] 


2 


2 


[1,2,1,2,1] 


2 


5 


3 


[1,2,2,1] 


2 


2 


[1,2,1,2] 






4 


[1,2,1] 


1 


2 








5 


[1,1] 


2 


2 








6 


[1] 


1 


1 








7 


□ 
















= 1 




v(u) = 


2/3 





In the following example, Lemmas 15. 1 1 and [5T2l are exploited to explicitly compute E[A^]. 

Example 4: Let us consider p = 4. The total number of partitions of V = {1, 2, 3, 4} is equal 
to .6(4) = 15. Considering the tree of depth p, we apply to each path the rules of Lemma [5T2l 
and we find that 14 paths (partitions) out of 15 reduce to the empty path, thus contributing 
with v{u) = 1. The only path that cannot be further reduced is u = [1, 2, 1, 2]. Thus, applying 
Lemma [54] with n = 2, we obtain v(u) = 2/3. From ((24b and considering all contributions, 
we obtain: 

E[Xy = (3 3 + (6 + (2/3) d ) /? 2 + 6/? + 1. 



VI. Convergence to the Marcenko-Pastur distribution 

In Section [V] we have shown that the moments of the asymptotic eigenvalues of are polynomials 
in (3, given by d24l) . In particular, the p-th moment E[A^ J has degree p — 1 and is given by the sum of 
B(p) positive contributions of the form v{u) d f3 p ~ k ^ . Since < v(u) < 1 and j3 > 0, for any d, the 
following inequality holds: 

i.e., for any given p and (3, the moments of the asymptotic eigenvalues decrease as the field dimension 
increases. The series E[A^ J, as a function of d, is positive and monotonically decreasing, thus it converges 
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to: 



nKo,p}= Km E[X p d J (27) 
Lemma 6.1: The moments E[A^ a] are the Narayana polynomials, given by 



k=l 

where T(p, k) = ^ Q! Z\) {k-i) are the Narayana numbers [34], [35]. Moreover, the random variable ^ 
follows the Marcenko-Pastur distribution [6] with pdf (see Figure [2): 



where ci, c 2 = (1 ± \//3) 2 , < /3 < 1, c 2 < x < c\. 
Proof: The proof is given in Appendix |D] 




Fig. 2. Marcenko-Pastur distribution 

In the following, we apply our findings to the study of the LMMSE of a reconstructed multidimensional 
field; in particular, we exploit the Marcenko-Pastur distribution to compute the expectation in (ITTb . 
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VII. Study of the reconstruction quality through the Marcenko-Pastur distribution 
Recall that the MSE provided by the LMMSE filter is [5]: 



MSE 



LMMSE 



E 



where X^g is distributed as the asymptotic eigenvalues of T^, with pdf fd,p(x). 
By using the Marcenko-Pastur distribution foa,f3 instead of f^p, we have: 



MSE r 



E 


a(3 


_ f Cl aPfooA*) 




J C2 x + af5 


-e + ^/e 2 -4f3 






2(3 





dx 



(30) 

zp 

where = 1+0(1 + a). 

Equation (l30l provides an approximation to the MSEqo, which, as shown in the following plots, can 
be exploited to derive the quality of the reconstructed field, given a finite d. 

We first consider d = 1 and compare in Figure [3] the expression of the MSEqo as in (fTTT) (solid 
lines) with the one obtained by using the Marcenko-Pastur distribution (dashed lines). The results are 
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presented as functions of the SNR m and for different values of j3. We computed (TTTT > by averaging over 
the eigenvalues of 200 realizations of the matrix Tx, with M = 150. The plot shows that, for small 
values of j3, the Marcenko-Pastur distribution (l30b yields an excellent approximation to the MSEqq, 
already for d = 1. Instead, for values of (3 greater than 0.2, the expression in 001 fails to provide a valid 
approximation. 

However, it is interesting to notice that, for d > 1, it is possible to obtain an accurate approximation of 
the MSEqo using the Marcenko-Pastur distribution, even for large values of (3. This is shown by Figures 
0] and [51 which plot the results obtained through (fTTT) and (l30l for (3 equal to 0.4 and 0.8, respectively. 
The results are presented as the SNR m varies and for different values of the field dimension d. 

Looking at Figure [4] we note that our approximation is tight for d > 2, while Figure [5J shows that, 
when d = 3, we still get a fairly good approximation for (3 as large as 0.8. 



April 22, 2008 



DRAFT 



20 




VIII. Conclusions 

We considered a large-scale wireless sensor network sampling a multidimensional field, and we 
investigated the mean square error (MSE) of the signal reconstructed at the sink node. We noticed that 
an analytical study of the quality of the reconstructed field could be carried out by using the eigenvalue 
distribution of the matrix representing the sampling system. Since such a distribution is unknown, we 
first derived a closed-form expression of the distribution moments. By using this expression, we were 
able to show that the eigenvalue distribution of the reconstruction matrix tends to the Marcenko-Pastur 
distribution as the field dimension tends to infinity. We applied our results to the study of the MSE of the 
reconstructed field, when linear filtering is used at the sink node. We found that, by using the Marcenko- 
Pastur distribution instead of the actual eigenvalue distribution, we obtain a close approximation to the 
MSE of the reconstructed signal, which holds for field dimensions d > 2. 

We believe that our work is the basis for an analytical study of various aspects concerning the 
reconstruction quality of multidimensional sensor fields, and, more generally, of irregularly sampled 
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signals. 

Appendix A 
The constraints 

Let us consider a vector of integers q of size p partitioning the set V = {1, . . . ,p} in k subsets Vj, 
1 < j < k, and the set of k constraints (|2TI ). We first show that one of these constraints is always 
redundant. 



A. Redundant constraint 

Choose an integer j, 1 < j < k. Summing up all constraints except for the j-th, we get: 

k k 

c h = Ch ~ c i = J2 £i - £ li+n - c i = ~ c i 

h=l h=l ieV 

which gives the j-th constraint since J2iev ^* — £ U+l] = ^- Thus, one of the constraints in (1211) is always 
redundant. Next, we show that the remaining k — 1 constraints are linearly independent. 

B. Linear independence 

The k constraints in (|2TT ) can be arranged in the form: ~W£ = with £ = [£i, . . . , £ p ] T and W being 
a k x p matrix defined as 

W = W / -W" (31) 



where 



(W) 



+1 i G Vj 



otherwise 

and W" is obtained from W' by circularly shifting the rows to the right by one position. Since one of 
the constraints (|2TT) is redundant, the rank of W is: p(W) < k — 1. Now we prove that the rank of W 
is equal to k — 1. 

Since the subsets Vj have empty intersection, the rows of W are linearly independent; hence, W 
has rank k. Also, W" is obtained from W by circularly shifting the rows by one position to the right, 
thus W" can be written as W" = W'S where S is the p x p right-shift matrix [36], i.e., the entries 
of the i-th row of S are zeros except for an entry equal to 1 at position [i + 1]. As a consequence, 
W = W — W'S = W'(I P — S) where the rows of the matrix I p — S are obtained by circularly shifting 
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the vector [+1, —1,0, ... ,0] and thus has rank p(I p — S) = p — 1. Hence, using the properties of the 
rank of matrix products reported in [36], we have 

p(W) = p(W'(I p - S)) > p(W') + p(I p -S) — p = k — l 

We recall that the system of linear equations Wz has a finite number of integer solutions bounded in 
[— M, . . . , M] d . The number of solutions decreases as p(W) increases. 

Appendix B 
Proof of Lemma IsTTI 

Proof: Using ([20]> and ([25]>, we obtain: 



fe(w) 

v{u) = lim — —. r— > 1 I 5(cj) 

v ^ 7 ^e£i i=i 

We first notice that nj=^ ^( c j) = S(W£) where the k(u>) x p matrix W is defined in Appendix IA1 
and S(W£) is a multidimensional Kronecker delta [33]. Since the rank of W is p(W) = k(u) — 1 (see 
Appendix lAl. then 5(W£) defines a subspace of IP with p—k(u)) + l dimensions. Therefore, considering 
that £ is a vector of integers with entries ranging in the interval [— M, . . . , M] and taking the limit for 
M — ► oo, we obtain u(u;) = Jj_i/2 1/2)* (Wz) dz where zeF and the function <5^ represents the 
Dirac delta. We have that ^(Wz) can be factorized as 

fc(w) 

^d(Wz) = J] <J d (rJz) (32) 
i=i 

where rj is the j-th row of W. As already shown in Appendix |A] one of the constraints (|2TT > is redundant 
and, hence, one of the factors in the right hand side of (l32l) . say the n-th, must not be included in the 
product. Now, moving to the Fourier transform domain, we can write: 5d{rjz) = exp(j27ryrjz) dy. 
Therefore, 

r fc(«) 
v(u) = / [] ^(rjz)dz 



1 [ / e^^d^dz 

,-1/2,1/2)* ^^./-oo 

/ / e>' 27r * r ? z dy„ dz 

V[-i/2,i/2)p Jm k - 1 
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where y n = [y u . . . , y n _i, y n+1 , y k(t<JJ )] T . Integrating first with respect to z, we get 

v{u) = f f e> 2 ^*=^ y " w *dzdy n 

jR fc -! J [-1/2,1/2)' 

V 



[ f f[ei 2 ^>'dzdy n 

il" J[-1/2,1/2)p ~^ 
f P /-1/2 

/ 11/ ei 2 ^>^d^dy n 

Jl'- 1 " 7-1/2 

J]^> 4 sinc(yT W4 )dy n 

i=i 

r p 
/ e i^^- w 'r[sinc(yTw i )dy n 



where w, is the i-th column of W, taken after removing its n-th row. By definition, the j-th rows of W 
and of W" contain both p- \Tj\ "0" and \Tj\ Since W = W - W", we have YH=i w * = and 

v(u)) = J Rfc _ 1 n?=i sin c(yn w dy n - Notice that, by definition of W (see <ED), y^w; = yj - y r \ Vn=Q 
if i 6 and [i + 1] € P,v. Moreover, by the definition in ( fl6l ), we have j/j = y Wi when i £ Vj. Thus, 



v ( u ) = / TT sinc (^ L=o dy n 

i=1 



Appendix C 
Proof of Lemma [5T21 (Simplification Rules) 

Let Vj be a singleton with Vj = i and u>i = j. We first notice that, since Vj is a singleton, [i— 1], ^ 
Vj. By applying Lemma I5TT1 with an arbitrary n ^ j, we have 

f P 

= / FT sinc(y Wh - | Vn =o dy n 

f P 



•sinc(y Wi - y W[<+1] ) | y „= dy n 
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We now integrate with respect to yj, with j = U{ and we obtain 



, p 
>(«)=/ Ff sinc(y, 

JR k - 2 



h=l 
h*[i-l] 



■ / sinc(y a , [ ._ 1] - yj)smc(yj - y Uli+1] ) \y n=0 dyjdy' n 
jm, 

f P 



h=l 
h^[i-l] 



■sincfe^, L=o dy^ 

,. p-i 

/ II sinc (^ - y^L+J l»»=o dy' n = v{u') 



where y' n and u>' have been obtained from y n and a; by removing their j-th and i-th element, respectively. 
Obviously y' n has size k— 1 and u/ has size p— 1. Let Vj be such that: "Pj = i, [i + 1], i.e., = = j. 

Then, 



V ( U )= IT sinC (^ - yw[fc+u) L=° d yr. 

= / FT sinc(y Wh - y u ) 



•sinc(y Wj | yn =o dy n 

^ [h+1] )sinc(y j - y^) | Wb=0 dy n 



/ TT sinc(?/ a 



/ IT sinc (^ - y^i h +i]) L=o dy r , 



h=l 
p-1 



/ JT sinc(y^ - yu') |j/„=o dy n = v(J) 
JK h=i 



where u' has been obtained from u by removing its z'-th element. 

Appendix D 
Proof of Lemma I6T1 

In order to prove Lemma 16.11 we first note that U Pt k may contain both crossing and non-crossing 
partitions [37]. 
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a) Non-crossing partitions: Every non-crossing partition contains at least a singleton or a subset 
with adjacencies, and therefore can be reduced by using the rules in Lemma 15.21 After reduction, the 
resulting partition is still non-crossing, thus it can be further reduced till the empty set is reached. It 



b) Crossing partitions: Recall that, in general, the coefficient v(u>) defined in (l26l) can be obtained 
by counting the solutions of the system of equations: W£ = where the fc(u;) x p matrix W contains 
the coefficients of the k(uj) constraints in (f2TT) . 
If cj 6 Qpfi is a crossing partition, then 

• k(u) > 2 (by definition, a partition with fc(aj) = 1 is always non-crossing) 

• it contains at least two subsets Vj and Vy, with j ^ f, which are crossing. 

Some crossing partitions can be reduced by applying the rules in Lemma 15.21 but, even after reduction, 
they remain crossing. 

Let us now focus on the crossing subset Vj of a partition u> which has been reduced by applying the 
rules in Lemma |5T2l Without loss of generality, we assume that \Vj\ = h, i.e. the partition Vj contains h 
elements with h > 2 since Vj is not a singleton. Then, by definition of the matrix W (see Appendix |A]) 
its j-th row, rj, contains h entries with value 1, h entries with value —1 and p — 2h zeros. We then 
build the 2 x p matrix W as W = [rj, — rj] T . Notice that W has rank 1 and the system of equations 
W£ = contains the constraints induced by a partition u> = [1, 2, . . . , 1, 2] with 2h entries. Since the 
system of equations W£ = contains a reduced set of constraints with respect to W£ = and, thus, a 
larger number of solutions, it follows that v{u>) < v{u>). 

It is straightforward to show that for a partition such as u>, with k(u) = 2, the coefficient v(u>) is 
given by Lemma [BTTI as v(u>) = L sinc(y) 2?i dy. This is a decreasing function of h and since h > 2 we 
have: 



follows that the non-crossing partition cj € Pi fc contributes to the expression of E[A^ J with a coefficient 
= 1. 




Therefore, we conclude that v(u) < v(u) < 1. 
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c) Crossing and non-crossing partitions: Let J)£ k , fi™ k C be, respectively, the set of crossing 



and non-crossing partitions of P) fe, with Vt c p k n 0™ fe = and n^UO^ = £l Pj k- Then, 



fc=i ujeQp.k 




fc=i 



where the equality (a) is due to the fact that for non-crossing partitions v(u) = 1, while for crossing 
partitions v(u>) < 1 and, hence, limd_> +00 v(u>) d = 0. In [38], it can be found that the number of 
non-crossing partitions of size A: in a p-element set is given by the Narayana numbers T(p, k) = |0™ fc | 
and therefore EfA^^] = Y%=iT{p,k)[3 p ~ k are the Narayana polynomials. In [6], it is shown that the 
Narayana polynomials are the moments of the Marcenko-Pastur distribution. 

Appendix E 
Proof of (fTTb 

We show that when the LMMSE filter is used, the expression of the asymptotic MSE is given by (fTTb . 
Indeed, by using (fTOt . (T7]), (O, and © we have: 



where = R^ + al and = G^G^. Substituting dU in the above expression and assuming E[aa^] = 
cr^I and E[nn^] = cr^I, we get 



where = (3Hd- Let us consider an analytic function g(-) in M + . Let X = UAU^ be a random positive 
definite Hermitian n x n matrix, where U is the eigenvectors matrix of X and A is a diagonal matrix 
containing the eigenvalues of X. By using the result for symmetric matrices in [39, Ch. 6] combined 
with the result in [40, pag. 481], we have: lim n ^oo iTrE[<?(X)] = E[g(A)] where the random variable 




T E [HA^Gdp-all 2 ] 

= Tr {(A^Rrf - IXA-iRd - I)t + aA^R^A, 



-i 



Tr {a(R d + al)- 1 } = Tr {a/3(T d + a/31)- 1 } 



A 
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A is distributed as the asymptotic eigenvalues of X. It follows that 



TrEhSCl^ + a/?!)- 1 ] 



a(5 




(2M + l) d 
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